Phase transition to chaos in complex ecosystems with non-reciprocal species-resource interactions

Non-reciprocal interactions between microscopic constituents can profoundly shape the large-scale properties of complex systems. Here, we investigate the effects of non-reciprocity in the context of theoretical ecology by analyzing a generalization of MacArthur’s consumer-resource model with asymmetric interactions between species and resources. Using a mixture of analytic cavity calculations and numerical simulations, we show that such ecosystems generically undergo a phase transition to chaotic dynamics as the amount of non-reciprocity is increased. We analytically construct the phase diagram for this model and show that the emergence of chaos is controlled by a single quantity: the ratio of surviving species to surviving resources. We also numerically calculate the Lyapunov exponents in the chaotic phase and carefully analyze finite-size effects. Our findings show how non-reciprocal interactions can give rise to complex and unpredictable dynamical behaviors even in the simplest ecological consumer-resource models.


D. Simulation and numerical methods
1. Numerical integration 2. Surviving/depletion cutoffs, numerical instability, and small immigration 3. Numerically solving the cavity self-consistency equations and calculating phase boundaries The objective of the cavity calculation is to find the steady-state behavior of the aMCRM.In particular, we will find the distribution of the steady-state abundances of the species and resources.The cavity method takes advantage of the system's self-averaging behavior.A system is said to have self-averaging behavior if the distribution of properties of constituents are independent of the exact realization of quenched disorder.This means that constituents' properties can be treated as random variates drawn from a distribution that is common to all systems with the same distribution of quenched disorder.A consequence of this is that an average of some observable taken over constituents given a fixed realization of quenched disorder is equal to the average of the observable for one constituent taken over all realizations of quenched disorder.The emergence of this self-averaging behavior is a consequence of the central limit theorem and quenched disorder, as we will see.In Fig. A5, we compare the distribution of steady-state resource and species abundance for a single realization of quenched disorder to the distribution of abundances predicted by the cavity method.

Setup
We begin by introducing the constituent averages, where we have introduced the notation X to denote the steady-state value of a quantity X.With these quantities defined, we can write the steady-state aMCRM as, where (A4)

Cavity solution
The cavity method begins by introducing a new species i = 0 and a new resource α = 0. We will later treat the new terms introduced by these variables as small perturbations.With the new species and resource, we can write the steady-state aMCRM for the existing species i = 1, . . ., S and resources α = 1, . . ., M as, For the new species and resource, the steady-state aMCRM says, Next, we will analyze the perturbed system relative to the unperturbed system.A quantity with \0 represents the value before adding the new variables to the system.Looking to Eqs.A5 and A6, we see that the presence of the new species and resource effectively perturbs the model parameters as, In the thermodynamic limit where M and S are large, we model the perturbation using linear response: where we have defined the susceptibility matrices: ) a. Self-consistency equations for species populations With the perturbation introduced and susceptibilities defined, we will exploit the linear response approximation to derive the self-consistency equations for the species populations.Substituting Eqs.A10, A11 into the aMCRM where we have defined q . Due to the self-averaging nature of the system and that the size of the perturbation is of order O(M −1/2 ), we have q R ≈ ⟨R 2 0 ⟩.Let Z N ∼ N (0, 1) be a unit normal random variable.The large-M limit approximate steady-state condition for the perturbing species becomes, 0 = N 0 g − σ c σ e ρχN 0 + σ g Z N . (A20) Additionally, observe that if χ → 0, the cavity solution diverges; this will be discussed in further detail in SI section A 5.
b. Self-consistency equations for resource abundances Now, we repeat this process to find a self-consistency equation for the resources.We substitute the linear response approximation for species into the aMCRM steady-state equation for the additional resource: Observe that the fourth term (involving χ iβ ) has zero mean and variance of order O(M −1 ); this can be seen by recalling that d i0 , d 0β , x i0 , x 0β are all independent for i, β ≥ 1.Similarly, we see that the variance of the fifth term (involving ν ).We will ignore fluctuations of order O(M −1/2 ).The mean of the fifth term is, where we have defined, which is the trace of the N i ← m j susceptibility matrix divided by S. Discarding terms of order O(M −1/2 ), we obtain, Now, observe the third and last terms are a sum of many independent random variables, meaning we can apply the central limit theorem and model the sum of these terms as a normal random variable.The mean of these terms is, The variance is, Var where q N ≡ 1 S S i=1 N 2 i\0 which is approximately equal to ⟨N 2 0 ⟩.The approximate steady-state condition for the added resource then becomes, where Z R ∼ N (0, 1) is a standard normal random variable.Solving for R 0 and discarding invadable solutions gives,

ReLU function-transformed normal distributions
In these computations, we regularly work with normal distributions that are transformed by the 'ReLU' function: where, is the standard normal CDF.The jth (j ≥ 1) moment is then, where 1 F 1 is the confluent hypergeometric function of the first kind, and Γ is the gamma function.Observe that W j (µ/α, σ/α) = α −j W j (µ, σ).Additionally, It follows from integration by parts, so the jth moment (j ≥ 1) is, ) 4. Final self-consistency equations Some essential quantities of interest are the expected fraction of surviving species ϕ N and fraction of non-depleted resources ϕ R .These quantities are computed using the moments calculated in section A 3 and Eqs.A20, A28: where ∆ g = g/σ g and ∆ κ = κ/σ κ and Θ is the Heaviside step function with the convention Θ(0) = 0. Next, we can differentiate our expressions for N 0 and R 0 to find, We can solve these two equations for χ, ν to obtain the relations, Next, we use Eqs.A20 and A28 and invoke our assumption that the system self-averages to find, The equations A39, A40, A41, A42, A44, A45, A46, A47 constitute the cavity self-consistency equations for the aMCRM model.They are eight independent nonlinear equations to solve for eight variables: ϕ N , ϕ R , χ, ν, ⟨N ⟩ , ⟨R⟩ , q N , q R .These equations are solved numerically using nonlinear least squares as discussed in SI section D 3.

Infeasibility of the cavity solution
Next, we will investigate when there may exist a solution to the cavity self-consistency equations.Observe that Eq. (A41) implies that as χ → 0, ν diverges; Eqs.A21, A29 indicate that N 0 and R 0 become singular.This is when, numerically, there fails to exist a solution to the cavity self-consistency equations.From Eq. (A42), we see that χ → 0 when, The cavity solution no longer has a solution when the ecosystem approaches the bound set by the principle of competitive exclusion.According to this calculation, ν diverging here means that if a species becomes marginally more fit when all niches are packed, it will disrupt the ecosystem.This relation defines a ninth equation for eight variables, so it determines a co-dimension-one boundary in the space of parameters.This boundary can be calculated numerically using nonlinear least squares as discussed in SI section D 3. When solving the cavity self-consistency equations numerically, we will see that there is a region of parameter space beyond this boundary where the least squares objectives are large compared to machine error (see Fig. D15 in SI section D 3).This boundary corresponds to the transition to the "unbounded growth" phase in the Lotka-Volterra literature (see Ref. [31]) because ⟨N ⟩ → ∞ as χ → 0 as can be seen in Eq. (A44).There is no unbounded growth in this model due to the negative feedback structure of consumer-resource models.In the cavity solution, this can be seen through ⟨R⟩ → 0 as ν → −∞ in Eq. (A45).That is, species abundances are predicted to diverge while abundances of all resources become zero.Additionally, as χ → 0, q R → 0, so σ g = σ m and g = −m, meaning ϕ CSI N = 1 − Φ(m/σ m ); therefore ϕ N and ϕ R approach nonzero values at this boundary.

Comparison of cavity and simulation results
For higher reciprocity levels ρ when the system is in the stable phase, the analytic solution produced from the cavity method matches the simulation results remarkably well.In Fig. A6, we compare the cavity results to numerical simulations for various 0 ≤ ρ ≤ 1 and fixed σ c = σ e = 3.5, corresponding to the gray dashed line in Fig. 3(a).For all predicted quantities, the chi-squared statistic is near the number of degrees of freedom, indicating that the cavity solution is a good fit to the simulation results.The cavity susceptibilities χ, ν for this slice are shown in Fig. A7.In this figure, we see that the infeasibility boundary occurs when χ → 0 and ν → −∞.
In Fig. A8, we compare the cavity results to numerical simulations for the various values ρ and σ c = σ e ≡ σ shown on the grid in Fig. 3(a).Again, the analytic solutions produced from the cavity method match the simulation results remarkably well in the stable phase.Cavity susceptibilities χ, ν for this grid are shown in Fig. A9.

Cavity susceptibilities
The cavity susceptibilities, have physical significance: χ and ν are the average linear-order response of a resource's abundance to a small change in its carrying capacity and a species' population to a small change in its natural mortality rate, respectively.These susceptibilities could be measured numerically by introducing small, nonzero-mean perturbations to K α and measuring M α=1 ∆R α / M α=1 ∆K α when the system is self-averaging.At steady state, the full susceptibility matrices can be computed using just the knowledge of which species and resources.Consider that the surviving species and resources satisfy the following matrix equation: where X * represents a matrix or vector with rows and/or columns corresponding to species and resource which are not surviving removed.This block-matrix equation can be solved to find, A9. Heatmap of cavity susceptibilities χ, ν as defined in lines A16, A24 for various values of σ and ρ, as in Fig. A8.
Given that perturbations are sufficiently small so that no species or resources go extinct and no new species or resources can invade, the susceptibility matrices are, [ν This means that when a steady state exists, the cavity susceptibilities can be computed exactly with knowledge of which species and resources are surviving.Notice that (e * ) T (c * (e * ) T ) −1 c * is an oblique projector because it is idempotent.The trace of this matrix is the dimension of the range of the projector.Therefore, the infeasibility boundary in the cavity solution occurs (or equivalently, the niches are fully packed) when the projector is full-rank.
Appendix B: Stability phase transition in the thermodynamic limit In nonlinear dynamics, a common approach to assess the appearance of instability is to begin with an assumption of stability and see what breaks down when the system is perturbed.In the case of the aMCRM, when there is stability, the results from the cavity solution will be valid, so we will use these results to determine the critical threshold for instability.
We consider perturbing all surviving species and resources, N are independent random variables all with mean zero and variance one [31].By considering perturbations on surviving species and resources and using cavity results (c.f., Eqs.A29, A21), we are assuming that we are working with an uninvadable steady state.From Eqs.A17, A25, for surviving species and non-depleted resources, Differentiating with respect to ε yields, dN dR are all independent, ⟨dR + 0 /dε⟩ = 0; similarly, ⟨dN + 0 /dε⟩ = 0.The first moment of the response to the perturbation does not give useful information about the stability; therefore, we turn to the second moment.First, we square the above quantities: Averaging over all sources of randomness and using , which follows from the assumption that the system self averages, we obtain: Heatmaps of the variances of the response of a surviving species and resources to a random perturbation (Eqs.B15, B16) as a function of the parameters σ and ρ with other parameters matching those in Fig. 3(a).The instability boundary is shown as a gray solid line, and the infeasibility boundary is shown as a gray dashed line. (B14) Solving this system of equations yields: These susceptibilities are the variance of the response of a surviving species or resource when the system is subject to a random perturbation.The divergence of these variances represent the breakdown of the mean-field approximation.Further, when these susceptibilities diverge, the uninvadable fixed point predicted by the mean-field approximation becomes dynamically unstable.Therefore, when these susceptibilities are divergent, whenever invasion is attempted, the system will exhibit dynamical fluctuations.We call the boundary in parameter space at which these susceptibilities diverge the instability boundary.These susceptibilities diverge when 0 = ρχ 1 − νρσ c σ e γ −1 2 − γ −1 ϕ N ϕ R .This condition along with the cavity self-consistency equations (Eqs.A39, A40, A41, A42, A44, A45, A46, A47) determines a co-dimension-one boundary in the space of model parameters.By using the relations in line A43, we obtain the following relation for the boundary at which the system becomes unstable: ⇓ instability boundary occurs when: (ρ ⋆ ) 2 = (# of surviving species) ⋆ (# of non-depleted resources) ⋆ , (B18) where X ⋆ denotes the value of a quantity X at the instability boundary.We are able to determine the critical level of asymmetry, ρ ⋆ , at which the ecosystem becomes unstable by solving a nonlinear least squares problem (see SI section D 3).The variances of the susceptibilities ⟨ dN + 0 /dε 2 ⟩, ⟨ dR + 0 /dε 2 ⟩ are shown for various values of σ = σ c = σ e and ρ in Fig. B10 with the instability and infeasibility boundaries overlaid.These variances of susceptibilities additionally have an asymptotic power law dependence on the distance from the instability boundary, |ρ − ρ ⋆ |, as shown in Fig. B11.The transition to the dynamic phase coincides in method to the transition to the "multiple attractors" phase discussed in the random generalized Lotka-Volterra model literature.This was first described in Ref. [31].
Appendix C: Chaotic nature of dynamics In the dynamic phase of the aMCRM, when an ecosystem is sufficiently large, the dynamics are chaotic.In this section, we provide numerical evidence that the dynamics of the aMCRM are chaotic in the classical sense (i.e., sensitive dependence on initial conditions) and in the sense of unpredictability (i.e., no clear patterns in the dynamics).As shown in Fig. 4(b) in the main text, the trajectories of the abundances of the surviving species and resources diverge from one another, indicating that the dynamics are chaotic in the classical sense.Evidence of chaos and unpredictability in the dynamic phase is present in other numerical experiments, as well.
In Fig. C12, we show that the dynamics of the aMCRM in the dynamic phase are unpredictable in the sense that a trajectory does not display any clear patterns in the projection onto the first three principal components of the correlation matrix of the time series of the abundances of the surviving species and resources.Such a projection is necessary to visualize the dynamics of the system, as they are very high-dimensional.

FIG. C12.
Projection of M + S = 256 + 256 = 512-dimensional dynamics onto the 3 highest-ranked eigenvectors of the correlation matrix (first 3 principal components) of the time series of the abundances.The lack of clustering of the trajectories in the projection indicates that the dynamics are chaotic and unpredictable.The parameters correspond to the gray star in Fig. 3(a).Before plotting, the simulation was run for 10 4 time units to eliminate any potential transients.

Analysis of the Lyapunov exponents
In Fig. 4(a), we show the maximal Lyapunov exponent for simulations classified by whether they reach a steady state for various values of ρ.To determine the Lyapunov spectrum, we use the lyapunovspectrum function in the ChaosTools.jlJulia package which employs the 'H2' method of Geist, originally stated in .This algorithm is described in detail in SI section A of Ref. [41], and its applications are described in Chapter 3.2.Conceptually, in order to compute the kth largest Lyapunov exponent, the algorithm evolves n ≥ k deviation vectors in the tangent space of the system and computes how the shape of the n-dimensional parallelepiped spanned by the deviation vectors evolve; the eigenvalues of the matrix describing the evolution of the parallelepiped are asymptotically related to the Lyapunov exponents.The Jacobian of the system when the dynamical variable is considered as (N 1 , . . ., N S , R 1 , . . ., R M ) is: The Lyapunov spectrum is computed with M = S = 512 species and resources, and the initial deviation vectors are chosen to be 8 unit vectors in the direction of randomly chosen species and resources.The parameters passed to the lyapunovspectrum function are the number of steps, N = 1000, the time step, ∆t = 5, and the time to wait before starting to record the Lyapunov spectrum, Ttr = 1000.The initial conditions are N i (0) = 1/S and R α (0) = 1/M .As seen in Fig. 4(a), the maximal Lyapunov exponent for a simulation that does not reach a steady state is positive while the maximal Lyapunov exponent for a simulation that does reach a steady state is negative.A simulation is categorized as reaching a steady state if the mean absolute value of all derivatives of the abundances at the end of the simulation is less than 10 −6 .These results indicates that the dynamics are chaotic in the classical sense for the simulations which have persistent fluctuations.The impact of system size on the probability of reaching a steady state is discussed in SI section E.

Analysis of the generalized alignment index (GALI)
FIG. C13.Generalized alignment indices (GALIs) for simulations in the dynamic and stable phases with M = S = 512 species and resources.The simulations all have the same sampled parameters; only ρ and order of the GALI are varied.For the simulations in the dynamic phase (solid lines), the GALI asymptotically decays exponentially with time, indicating chaos, while for the simulations in the stable phase (dashed lines), the GALI asymptotically approaches a nonzero value, indicating the dynamics achieve a steady state.Simulation parameters for the GALI in the stable and dynamic phases correspond to the gray dot and star in Fig. 3(a), respectively.
An alternative method to determine whether a system is chaotic is to use the generalized alignment index (GALI) [41][42][43].The GALI is a measure of how vectors in the tangent space of a trajectory align with each other.Let ŵ1 (0), . . ., ŵk (0) ∈ R M +S be linearly-independent unit deviation vectors that evolve according to, where x(t) = (N 1 (t), . . ., N S (t), R 1 (t), . . ., R M (t)) and J(x(t)) is the Jacobian matrix of the system in state x(t).These deviation vectors are normalized to have unit length at regular (small) time intervals.The order-k GALI is defined as, which is the volume of the k-dimensional parallelepiped spanned by the time-evolved deviation vectors.If the system is chaotic, the GALI will decay exponentially with time, indicating that the deviation vectors are becoming more aligned with each other due to the exponential divergence of nearby trajectories.If the system is not chaotic and reaches a steady state, the GALI will asymptotically approach a nonzero value, indicating that the deviation vectors form a nonzero volume in the tangent space because no single direction dominates the dynamics.Alternatively, if the system is asymptotically periodic or quasi-periodic and does not reach a steady state, the GALI will decay as a power law with time.
In Fig. C13, we show the GALI for simulations in the dynamic and stable phases for k = 2, . . ., 7 with M = S = 512 species and resources.We use the gali function from the ChaosTools.jlJulia package [42].The simulations all have the same sampled parameters; only ρ and the order of the GALI are varied.For the simulations in the dynamic phase (solid lines), the GALI asymptotically decays exponentially with time, indicating chaos, while for the simulations in the stable phase (dashed lines), the GALI asymptotically approaches a nonzero value, indicating the dynamics achieve a steady state.The initial deviation vectors are chosen to be k unit vectors in the direction of randomly chosen species and resources.The arguments passed to the gali function are the run time, T = 2e4, the time step, ∆t = 5., the threshold at which to stop the simulation, threshold = 1e-22, and the order of the GALI, k.

Analysis of effective dimension of dynamics
Chaotic systems can often be characterized by their effective dimension, which is the number of degrees of freedom that are relevant to the dynamics.One measure of the effective dimension is the Kaplan-Yorke (KY) dimension [44,57,58], which is the number of Lyapunov exponents that are positive.The KY dimension is the linearly interpolated point at which the cumulative sum of the Lyapunov exponents crosses zero: We estimate the KY dimension of the aMCRM using the Lyapunov spectrum measured using the methods described in SI section C 1 and the kaplanyorke dim function from the FractalDimensions.jlJulia package [58].We find that the KY dimension of the aMCRM in the dynamic phase for the parameters corresponding to the gray star in Fig. 3(a) is D KY = 15 ± 7, which is an average over 64 simulations.The number of surviving species and resources in these simulations at the end of the simulation is S ⋆ + M ⋆ = 74 ± 17.The ratio of the KY dimension to the number of surviving species and resources is D KY /(S ⋆ + M ⋆ ) = 0.21 ± 0.09.The number of species and resources that are transitioning between high-and low-abundance states is M t + S t = 38 ± 19, and the ratio of the KY dimension to the number of species and resources that are transitioning between high-and low-abundance states is D KY /(M t + S t ) = 2.6 ± 1.2 ≈ O(1).From this, we hypothesize that the effective dimension of the dynamics is approximately this number of species and resources 'jumping' between high-and low-abundance states.The intuition behind this phenomenon is that the hypothesis that dynamics of the aMCRM in the dynamic phase are dominated by the species and resources that are transitioning between high-and low-abundance states [26].We hope to explore this connection further in future work.Similarly, to find the instability boundary, we add a term to the objective function corresponding to the instability condition (Eq.(B18)): Minimizing these objective functions while allowing one additional model parameter to vary gives the infeasibility and instability boundaries shown in all figures in this paper.
In order to understand whether the aMCRM achieves a steady state in the thermodynamic limit, we must understand how finite system size and a finite simulation time impact simulation results.To assess the existence and dynamical nature of a steady state in a simulation, we numerically calculate a time to steady state (TtSS).After running a simulation, we iterate through the time series and find the time at which the mean absolute value of the derivatives last dips below a threshold, 10 −10 .The tolerance of the solver (Tsit5 [33,35]) is 10 −14 , so this threshold is well above the numerical error of the simulation.Because we simulate to a finite time T max = 2 16 + 10 4 = 75536, it is not possible to distinguish between a system that will reach a steady state at some time beyond T max and a system that will never reach a steady state.Furthermore, when performing these simulations, especially with large system sizes, errors occasionally occur in the solver, and the simulation is aborted; while we do not present finite-size scaling results in regions where the simulations are frequently aborted, a robust analysis should account for this.Here, we develop a custom maximum likelihood estimation (MLE) technique to fit the distributions TtSS's at different system sizes, taking into account these two situations.We then examine the scaling of the resulting fit parameters to determine their behavior in the thermodynamic limit.

Modeling simulation outcomes
To construct a MLE model, we first must model the statistical process which determines each possible outcome of a simulation.We observe that simulation outcomes fall into three categories: (i) The simulation reaches steady state within the maximum simulation time T max .
(ii) The simulation does not reach steady state within T max .
(iii) The solver encounters an error and the simulation aborts.
First, we represent the probability of a simulation encountering an error as p err .Next, we consider whether a simulation could possibly reach a steady state if given enough simulation time.We represent the probability of a simulation reach steady state in the infinite time limit as ϕ SS .Finally, we define the probability density function p SS (t) of TtSS's for those simulations that can reach steady state when given enough time.
Using these definitions, we can now construct the probabilities for each of the cases defined above.Let T be a random variable representing the observed outcome of a simulation with possible values T ∈ [0, T max ] ∪ {noSS, error}, where T takes a finite value in the interval [0, T max ] for case (i), or the special values noSS or error in cases (ii) and (iii), respectively.
To model case (i), we note that it results in an outcome represented by a continuous numerical value, while cases (ii) and (iii) represent categorical outcomes.For ease of derivation, we also convert T to a categorical outcome in the prior case by artificially breaking our simulation interval into small discrete time steps of size ∆t.Later, we will take the limit ∆t → 0, removing the dependency of our model on this step size.The probability of a simulation completing successfully and reaching a steady state in the time interval t < T < t + ∆t with T < T max is approximately For case (ii), simulations do not encounter an error, but do not reach steady state within T max .We must consider two separate possibilities: a simulation may never reach steady state, even if T max → ∞, or would reach steady-state after further simulation, T > T max .In the first case, the probability is simply (1 − p err )(1 − ϕ SS ).If a simulation would eventually have reached steady-state if we had continued simulating, we must use our probability density function p SS (t).However, we do not know exactly when the simulation would have finished, so we must sum across all discrete time intervals greater than T max .The total probability for case (ii) is then where we have taken the limit ∆t → 0 and simplified, defining F SS (T ) =

Scaling of parameters with system size and discussion
The simulated data for different system sizes in the stable and dynamic phases can be visualized by analyzing the empirical cumulative distribution functions (CDFs) of the observed times to steady state called T SS above.In Fig. E16, we show the CDFs for the stable and dynamic phases for 4 ≤ M, S ≤ 1024 along with the CDFs for the fitted distributions.A clear difference is apparent between the stable and dynamic phases, and the fitted distributions match the simulated data well.A scaling collapse of the CDF based on the fitted distributions is shown in Fig. E17.
The fit parameters for the stable and dynamic phases are shown in Fig. E18.In the stable phase, the estimated probability of reaching steady state, ϕ SS remains approximately equal to 1 for all system sizes, and the estimated timescale τ asymptotically increases with system size as a power law.In the dynamic phase, ϕ SS decreases with system size towards 0, and τ asymptotically increases with system size also as a power law, but with a larger exponent.In both phases, α approaches 1 asymptotically from above.
We provide fits of ϕ SS as a function of system size M using curves of the form ϕ SS (M ) = ϕ max −∆(1−exp(−(M/ξ) κ )).For the dynamic phase ϕ max = 1.0, ∆ = 1.0, ξ = 278.5, and κ = 1.5; for the stable phase, ϕ max = 1.0, ∆ = 0.0, κ remains approximately equal to the value at which it is initialized in the optimization algorithm, and ξ approaches the maximum value set in the optimization algorithm.
We also provide fits of the timescale τ as a function of system size M using power laws of the form τ (M ) = bM a .Curves are fit for data M, S ≥ 2 9.5 .In the stable phase, a = 0.4 and b = 580; in the dynamic phase, a = 0.9 and b = 70.
This finite size scaling is run in the range 512 ≤ M, S ≤ 1024 but is not shown in Fig. E18 because solver errors are occasionally present.More significantly, as very few simulations reach steady state in the dynamic phase when the system size is sufficiently large, the MLE procedure cannot clearly distinguish between cases where τ is very large or ϕ SS is very small; this ambiguity can be seen more clearly in Fig. E16.With the prior knowledge of the finite size scaling as shown in Fig. E18 where the MLE procedure succeeds and there are no solver errors, we conclude that at these larger system sizes, it is both the case that τ grows large (in polynomial order of system size) and ϕ SS approaches either zero or one.That fact both the timescale τ diverges and ϕ SS approaches either zero at large system size leads us to conclude that in the thermodynamic limit, all systems in the stable phase eventually reach steady state for long enough times, while in the dynamic phase, no systems reach steady state.Note that the curves collapse onto each other in both phases for system sizes that are moderately large M, S ⪆ 2 4.5 .This indicates that we have successfully identified the asymptotic scaling of the form of the TtSS distribution with system size.FIG.E18.Impact of finite size effects on the dynamics of the aMCRM.Rows correspond to different parameters describing dynamics: the estimated probability a simulation of a given system size will reach steady state (top) and the timescale of the time to steady state of simulations of a given system size that do reach steady state (bottom).Columns correspond to simulation data from different phases: the stable phase (left) and the dynamic phase (right); the parameters corresponding these simulations are marked in the phase diagram in Fig. 3 with a circle and star, respectively.The probability of reaching steady state in the stable phase remains at 1 for all system sizes, while the probability of reaching steady state in the unstable phase decreases to zero with increasing system size, indicating that stable dynamics are not observed within the dynamic phase in the thermodynamic limit.The timescale of the time to steady state increases with system size in both the stable and dynamic phases.The uninvadable steady state is found using an iterative constrained optimization method described in [16].Cavity parameters are found by solving the self-consistency equations using non-linear least squares as described in appendix D 3. All parameters are identical to those in Fig. 3(a) except that the random matrices and vectors are sampled as: This choice of distributions ensures that d iα , x iα , δK α , δm i are still standard random variables and appropriate thermodynamic scaling are used to obtain cavity results that are identical to those in Fig. 3(a).

E
example dynamics) simulation parameters Figure 3 (phase diagram) simulation parameters Figure 4 (Lyapunov dot plot and diverging trajectories) simulations details FIG. A5.Comparison of the empirical cumulative distribution functions (eCDFs) of the steady-state abundances of species (S = 512) and resources (M = 512) for a single realization of quenched disorder and the distribution predicted by the cavity calculation; simulations are performed in the stable phase.See appendix G for simulation parameters.
FIG. A6.Comparison of cavity results to numerical simulations for various ρ and fixed σc = σe = 3.5, corresponding to the gray dashed line in Fig. 3. (a) Average populations of species, ⟨N ⟩, and abundances of resources, ⟨R⟩.(b) Standard deviations of species populations, ⟨N 2 ⟩ − ⟨N ⟩ 2 , and resource abundances, ⟨R 2 ⟩ − ⟨R⟩ 2 .(c) Fractions of surviving species, ϕN , and non-depleted resources, ϕR.Points are means of numerical simulations, lines are cavity results.Error bars are standard deviations of numerical simulations; some error bars are smaller than the points.A dashed line is shown at ρ ⋆ = 0.72, the critical value of ρ at which the system transitions between the stable and dynamic phases.For each value of ρ, 64 simulations are performed with M, S = 512 resources and species.
FIG. A7.Cavity susceptibilities for the same parameters as in Fig. A6.The susceptibilities χ, ν as defined in lines A16, A24 are shown on the left and right, respectively.
FIG. A8.Heatmap comparison of cavity predictions and simulation results.Results from 32 simulations for various values of σ and ρ are displayed in the left column of each pane and the cavity predictions are displayed in the right column.The instability boundary is shown as a solid black line while the infeasibility boundary is shown as a dashed black line.For the cavity prediction plots, the parameter regions for which no cavity solution exists have a diamond pattern overlaid (c.f., Fig. D15 and SI section D 3).The parameters used are those used in Fig. 3(a).
FIG. D14.Log histogram of the species and resource abundances at steady state for a system with M, S = 256 and ρ = 0.9.There is a clear gap in abundance between species/resources that are surviving and those that are extinct.The immigration rate used is λ = 10 −10 .

T0
FIG.E16.Empirical cumulative distribution functions of T , the time to steady state, for simulations of different system sizes in the stable phase (left) and the dynamic phase (right).The maximum likelihood estimate distributions as described in SI section E for the time to steady state for the various system sizes are shown with dashed lines.The maximum value shown on the horizontal axis is the maximum time simulated, Tmax = 2 16 + 10 4 = 75536.
FIG. E17.Scaling collapse of TtSS CDFs.CDFs of observed TtSSs are scaled by the MLE-fitted value of ϕSS where the horizontal axis is scaled by the MLE-fitted values of τ , α.Note that the curves collapse onto each other in both phases for system sizes that are moderately large M, S ⪆ 2 4.5 .This indicates that we have successfully identified the asymptotic scaling of the form of the TtSS distribution with system size.

Figure F19 (
Figure F19 (phase diagram with uniform sampling) simulation parameters